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ABSTRACT 

Many physical properties of galaxies correlate with one another, and these correla- 
tions are often used to constrain galaxy formation models. Such correlations include 
the color-magnitude relation, the luminosity-size relation, the Fundamental Plane, etc. 
However, the transformation from observable (e.g. angular size, apparent brightness) 
to physical quantity (physical size, luminosity), is often distance-dependent. Noise in 
the distance estimate will lead to biased estimates of these correlations, thus com- 
promising the ability of photometric redshift surveys to constrain galaxy formation 
models. We describe two methods which can remove this bias. One is a generalization 
of the l^nax method, and the other is a maximum likelihood approach. We illustrate 
their effectiveness by studying the size-luminosity relation in a mock catalog, although 
both methods can be applied to other scaling relations as well. We show that if one 
simply uses photometric redshifts one obtains a biased relation; our methods correct 
for this bias and recover the true relation. 

Key words: methods: analytical, statistical - galaxies: formation — cosmology: 
observations. 



1 INTRODUCTION 

The 'configuration space' we use to describe galaxies is large, 
but galaxies do not fill it. The luminosity, color, size, surface 
brightness, stellar velocity dispersion, morphology, stellar 
meiss, star formation history and spectral energy distribu- 
tion of a galaxy are all correlated with one another. These 
correlations encode important information about galaxy for- 
mation, and so quantifying them provides important con- 
straints on models. 

Current (e.g. SDSS, Combo-17, MUSYC, Cosmos) and 
planned surveys (e.g. DES, LSST, SNAP) go considerably 
deeper in multicolor photometry than in spectroscopy, or 
are entirely photometric. For such surveys, reasonably ac- 
curate photometric redshift estimates are or will be made. 
The question then arises as to which galaxy observables and 
correlations are affected by the noisy distance estimate asso- 
ciated with a photometric rather than spectroscopic redshift. 
The most widely studied property is luminosity - clearly, er- 
rors in the distance result in incorrect luminosity estimates. 
If not accounted for, this leads to a biased estimate of the 
luminosity function (e.g. Subbarao et al. 1996). Hence, there 
has been considerable effort devoted to the question of how 
to correct for this bias (e.g. Chen et al. 2003), and the prob- 
lem has now been solved (Sheth 2007). 

The next step is to recover an unbiased estimate of 
not just the luminosity function, but the joint distribution 



of luminosity, color, size, etc., from photometric redshift 
datasets. The main goal of the present work is to provide 
an algorithm which does this for a magnitude limited pho- 
tometric redshift survey. Because the same distance error 
which leads to a mis-estimate of the luminosity will pro- 
duce a correlated mis-estimate of the size, we have chosen 
to phrase the discussion in terms of the size-luminosity re- 
lation - it exhibits all the features of interest. 

Section 2 illustrates the nature of the problem by show- 
ing the bias in the size-luminosity relation which results from 
treating photometric redshifts as though they were spectro- 
scopic redshifts. This is done by constructing a mock galaxy 
sample and then perturbing the true redshifts to mimic pho- 
tometric redshift errors. Section 3 places this problem in the 
more general context of inverse problems in statistical as- 
tronomy, and argues that a deconvolution algorithm, such 
as that due to Lucy (1974), is well-suited to removing the 
bias. It shows the result of applying this non-parametric de- 
convolution technique to a mock galaxy sample. Section 4 
provides a maximum-likelihood formulation and solution of 
the problem. A final section summarizes our findings and 
discusses possible further studies and applications. 

Where necessary, we write the Hubble constant as Ho — 
lOO/i km s~^ Mpc~^, and we assume a spatially fiat cosmo- 
logical model with {QM,^A,h) = (0.3,0.7,0.7), where 
and JIa are the present-day densities of matter and cosmo- 
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Figure 1. Distribution of spectroscopic and photometric red- 
shifts in our mock catalog which was set-up to mimic the SDSS 
early-type galaxy sample. Contours show levels which are 1/2" 
times the height of the maximum value, with n running from 1 to 
5. The photo-2 error distribution was assumed to follow equation 



logical constant scaled to the critical density. We use Dl (z) 
to denote the luminosity distance; the angular diameter dis- 
tance is Da{z) = Dl(z)/(1 -I- z)'^. 



2 CORRELATIONS WITH OBSERVABLES: 
EFFECT OF DISTANCE ERRORS 

In what follows, we use the luminosity-size relation to illus- 
trate how photo-z errors lead to biases. 

We begin by generating a magnitude-limited mock 
galaxy catalog with parameters chosen to mimic those 
of SDSS early-type galaxies in the g band, following the 
method given by Bernardi et al. (2003). The redshift range 
is restricted to the interval 0.01 ^ z ^ 0.3. We ignored pas- 
sive evolution of the luminosities and colors, as well as K- 
corrections. The simulated magnitude-limited catalog has a 
similar dN/dz distribution to that observed, and the distri- 
bution of apparent magnitudes, angular sizes, and velocity 
dispersions are very similar to those in the real data. 

We then model photometric redshifts by assuming 

that 



dx 



exp(— 7a;) 

r(7) ' 



(1) 



where x = Di,{(^)/ Dl{z) is the ratio of the photo-z based 
luminosity distance to the true one, and 7 = 5. This distri- 
bution has {x) = 1 and — I/7. With 7 = 5, this error 
distribution is substantially worse than typical photomet- 
ric redshift errors. Figure [1] compares ( and z. Note that 
the analysis which follows is not tied to this functional form 
for the photo- 2: error distribution; we are simply using it to 
illustrate our methods. 

In what follows, we use Al to denote the true abso- 
lute magnitude, and M to denote that estimated using 
rather than z. And, with some abuse of notation, we use 
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Figure 2. Distributions of intrinsic and estimated absolute mag- 
nitudes (top panel) and sizes (bottom panel) which result from 
the differences between true and photo-2 shown in the previous 
Figure. 



R to denote logj^p of the physical size 6Da{z), where 9 is 
the measured angular size and Da{z) is the angular diame- 
ter distance defined earlier. The estimated size based on the 
photometric redshift is 



7^ = log,o[eDA{0] = R + log,o[DA{0/DAiz)] 
= R-{M- M)/5 - 21ogio[(l + 0/(1 + z) 



(2) 



in principle, there are also evolution and K-correction terms 
which we set to zero. Figure [2] compares M with M and TZ 
with R. 

The qualitative nature of the distributions in the two 
panels are easy to understand. The distribution in M is 
broader than that in M, as is the distribution of TZ com- 
pared to R: photometric redshift errors have broadened both 
distributions. However, the changes to the estimated abso- 
lute magnitude and size are not independent. Assuming an 
object is closer than it really is makes one infer a fainter lu- 
minosity and smaller size than it really has. These correlated 
changes can have a dramatic effect on the size-luminosity re- 
lation, since photo-2: errors move each galaxy in the R — M 
plane left and down or right and up. In general, these mo- 
tions are not parallel to the principal axis of the true rela- 
tion, so the mean relation in the photo-z catalog, {7?.|A^), 
need not be the same as the true relation {R\M) . In our mock 
catalog, {TZ\M) oc -0.20 A^, whereas {R\M) cx -0.27 M 
(see Figure [3}. 
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Figure 3. Effect of photo-2 on the size-luminosity correlation in our mock catalog. In the left panel, contours and solid line show the 
TZ ~ M relation associated with photo-^s, whereas the right panel shows the intrinsic R — M relation. Note the strong bias (shallower 
slope in panel on left) which results from the fact that the photo-2; distance error moves points down and left or up and right on this 
plot. Squares in left panel show the binned starting guess for the 2d deconvolution algorithm, and triangles in right panel show the result 
after 20 iterations (equation (TSjl . Convergence to the correct solution is clearly seen. 



3 A NON-PARAMETRIC 

DECONVOLUTION-LIKE METHOD 

Sheth (2007) shows that to extend Schmidt's (1968) Vmax 
estimator of the luminosity function 4>{M) so that it pro- 
duces unbiased resuhs in photo-z surveys, it is helpful to 
think of J\f{M), the number of observed objects with esti- 
mated A4, as being a convolution of the true number which 
have M, N{M), with the probability that an object with 
magnitude M is thought to have magnitude M. The lumi- 
nosity function is then estimated by first deconvolving the 
distribution of J^{M) to obtain N{M), and then using the 
fact that N{M) = 4){M) Vra^^M). 

For the same reason, the present problem may be 
thought of as a two dimensional deconvolution problem (or 
an n-dimensional deconvolution if we were interested in the 
full manifold, rather than just two-dimensional projections 
of it). Specifically, let 

iV(M,i?) = N{M)p{R\M) = ^KnaxlM) (i,{M)p{R\M) (3) 

denote the (true) number of galaxies with absolute magni- 
tude M and size 7? in a magnitude-limited catalog. Here 
p{R\M) is the probability of having size R when the magni- 
tude is M. Similarly, set M{M, 7^) = M{M) p{n\M). Then 

N{M,1Z) = j AM j dR N{M,R) p{M,TZ\M,R) 

= j AM N[M,TZ~ F{M,M)]p{M\M), (4) 

where F = (X - M)/5 4- 21og[(l + C)/(l + z)\. Our prob- 
lem is to obtain a reliable estimate of the intrinsic N{M, R) 
given the photo-z biased J\f{A4,'R,) and the error distribu- 



tion p{M\A4). We do this using the deconvolution algorithm 
proposed by Lucy (1974). 

Before we present our algorithm and results, it is worth 
noting that we could have attempted to invert equation 
in other ways. Classical naive 'exact' inversion methods in- 
clude matrix-quadrature techniques (reduction of the inte- 
gral equation to a linear matrix system), polynomial ex- 
pansion methods, singular function expansion and product 
integration methods. These typically run into difficulties 
because the measured data function usually cannot pro- 
vide sufficient information on the high frequency compo- 
nents of the solution. A standard non-classical technique 
(Phillips 1962; Twomey 1963; Tikhonov 1963) is the method 
of regularisation. A "regularisation" parameter is intro- 
duced, which balances the size of the residual against the 
smoothness of the solution, and the problem is turned into 
one of minimization. However, there is no general strategy 
for choosing the optimum regularisation parameter; this led 
Lucy (1974) to formulate his algorithm, and is what has 
led us to choose his algorithm over these others. One might 
argue that we are likely to know a fair amount about the 
expected form of the intrinsic distribution (e.g., luminosity 
functions are rather smooth, and conditional distributions 
tend to be bell-shaped) , so it may be that these other meth- 
ods are worth pursuing further. This is the subject of work 
in progress. 



3.1 The deconvolution algorithm 

The general two-dimensional problem is that of estimating 
the frequency distribution ^'(^',??') of the quantities ^' and 
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r]' when the available measures Xi,y'i; x'2,y2', x'pf,y'pf are 
a finite sample drawn from an infinite population character- 
ized by: 



$(a;,y)= / d^ / d?? f?) p(a;, ??)• 



(5) 



Here ^{x, y) is the function accessible to measurement and 
pix,y\^,ri) is the conditional probability that x' will fall in 
[x, a; + da;] when it is known that ^' = ^, and that y' will fall 
in [y,y + dy] when rj' = rj. In many cases <1> and '3/ represent 
probability density functions, so they obey normalization 
and non-negativity constraints. 

The iterative procedure for generating estimates to 'I' 
presented in Lucy (1974) is 

V) = V) fdx [ dy p{x, v), (6) 



cE>'-(a;,y) 
where 

^"{x,y)= dU dr? *''(C,J7) 



(7) 



The index r indicates the rth iteration in the sequence of 
estimates, and $ is an approximation to $ obtained from 
the observed sample. Convergence is achieved if 'l?'^ = "!>. 
The starting approximation ^°(^,?7) should be a smooth, 
non-negative function having the same integrated density 
as the observed distribution. The extension to n-dimensions 
is obvious, although, if the number of dimensions is large, 
then performing the multi-dimensional integrations effi- 
ciently may become challenging. 

The two-dimensional problem simplifies if, as happens 
in our problem. 



P{x,y\^,v) 



p{x\^)5-D 



(8) 



because the delta function simplifies one of the integrals. 
The iterative scheme becomes 



■I'[x,77-F(x,0] 
cE>'-[a;,ry-F(s,e)] 



with 



p{x\0 (9) 



(10) 



Partial distributions can be easily computed via marginal- 
ization from 'I''""''^: 



/ dr,*'+i(e,7)). 



*'-+^(r,)= / de*'--^(e,77) 



r+l / 



and 



drj rj 



*'-+i(c) ■ 



(11) 



(12) 



(13) 



The generalization to n-dimensions is a straightforward ex- 
tension of the expressions above, so we do not present ex- 
phcit expressions. Note that, just as happens in the 2- 
dimensional case presented here, delta-functions will reduce 
the n-dimensional problem to a simple one-dimensional inte- 
gral, because the same distance error affects all n quantities. 
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Figure 4. Reconstruction of the intrinsic N{M) distribution from 
the distribution of estimated redshifts. Dashed histogram shows 
the observed absolute magnitude distribution, used as a starting 
guess. Jagged lines show the reconstructed intrinsic distribution 
after 20 iteration, using the simpler Id algorithm (blue) or a 2d 
iterative scheme (cyan). 
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Figure 5. Reconstruction of the intrinsic N{R) distribution from 
the distribution of estimated redshifts. Line styles same as previ- 
ous Figure. 



3.2 Results 

The formalism outlined above is readily applicable to the 
size- luminosity correlation if we interpret {x,y) as the esti- 
mated absolute magnitudes and sizes (A1,7^), and (^,7?) as 
the true intrinsic ones, {M,R). 

Figures |4] and [5] show how well this method recovers the 
intrinsic distribution of absolute magnitudes and sizes (solid 
histograms) . The broad dashed histograms show the photo-z 
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Figure 6. Examples of reconstructed conditional distributions p{R\M) bins in magnitude of width AM = 0.355 centred on AI = 
—22.046, —21.691, —20.982. Jagged lines show the distributions recovered by the 2d deconvolution algorithm after 20 iterations. 



derived A4 and TZ distributions. These were used as a con- 
venient starting guess in the deconvolution algorithm, al- 
though prior knowledge about the expected intrinsic shapes 
could have been used instead. The reconstruction after 20 
iterations is shown by the jagged lines. Of course, we could 
have chosen to reconstruct N{M) directly from Af{M), us- 
ing the Id deconvolution algorithm outlined in Sheth (2007) . 
This procedure converges in about 5 iterations to the distri- 
bution shown by the crosses. Note how well the deconvolved 
distributions match the intrinsic ones. 

A more stringent test is to check if the conditional dis- 
tributions, p{R\M), are also well recovered. Figure [S] shows 
p{R\M) for three bins of width AM = 0.355 centred at 
-22.046, -21.691 and -20.982. Clearly, the method works 
well. 

The means of these recovered distributions can be used 
as an estimate of the recovered size-luminosity relation: 
(R\M). Recall we had noted that this relation was rather 
strongly biased because use of the photo-z distance estimate 
means that the error in the size is correlated with that on the 
luminosity. The squares in the left hand panel of Figure [3] 
show the starting guess for our algorithm, and the triangles 
in the right hand panel show the reconstructed relation ob- 
tained from the deconvolution procedure — it is in excellent 
agreement with the true one. 



4 A MAXIMUM-LIKELIHOOD METHOD 

Sheth (2007) describes an algorithm which produces a 
maximum-likelihood estimate of the luminosity function 
from magnitude limited photo-z datasets. It is straightfor- 
ward to extend that analysis to the present case, in which 
the quantity of interest is not just the distribution of lumi- 
nosities, but the joint distribution of luminosity and other 
observed physical parameters. 

Let Adi denote the vector of physical quantities for 
galaxy i estimated using the photometric redshift when 



computing distances, and let J\f{M.i, (i\a) denote the num- 
ber of galaxies in a magnitude limited catalog that have es- 
timated redshifts (^i and estimated luminosities, sizes, etc., 
M-i, when the model for the intrinsic joint distribution of 
physical quantities is specified by the parameters a. Then, 
dropping the understood i-index dependence, 

J^{M.,C\a)= J dz^(l>{M\a)p((\z,m{M,(:)), (14) 

where p{(^\z,m) is the photo-z redshift error distribution, 
and cjf>(A^|a) denotes the true joint distribution of physical 
quantities, but evaluated at values which account for the 
fact that the same photo-z error which affects the absolute 
magnitudes affects the other observables. E.g., in the size- 
luminosity relation we have been considering. 



M = {M,'JZ) and M^{M,R), (15) 
where 

M = M~5\og,,[D^iz)/D^{0] (16) 

R = TZ + loe,^,[DA{z)/DA{0]. (17) 



The factor of m is the apparent magnitude associated with 
photo-z redshift and absolute magnitude A4; of course, m 
is the same for true redshift z and absolute magnitude M. 

The predicted number of objects with photometric red- 
shift ^ when the model has parameters a is 

Af{(\a) = J dM.Af{M,C\a). (18) 

If the redshift-error distribution is independent of m, then 

7V(C|a) = J dz{dV./dz)S{z,a)p{C\z) 

= J dzN{z\a)p{C\z), (19) 

where 
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/•M„i„{z) 

S{z,a) = / dM(^(M). (20) 

This shows that 7V((^|a) is just the convolution of the intrin- 
sic redshift distribution (in a flux-limited catalog) with the 
redshift-error distribution. 

The expressions above generalize those given in Sheth 
(2007), where M = M and M = M, with M given by 
equation (|16p . Hence, by analogy to when the distances are 
known accurately, the likelihood to be maximized is (rein- 
troducing the index i) 

L{a)=Y[p., where p. = r^^i^^L^. (21) 

i 

The analysis in Sheth (2007) can now be followed to show 
analytically that this is indeed the appropriate expression 
for the likelihood, so we do not reproduce it here. 



correlation between stellar velocity dispersion a and lumi- 
nosity L which arises if the distance indicator used to esti- 
mate L is correlated with cr (as may happen in the local Uni- 
verse, where peculiar velocities make spectroscopic redshifts 
unreliable distance estimators). It may be that the meth- 
ods presented here would allow an accurate reconstruction 
of the true relation from the biased one. This is the subject 
of on-going work. 
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5 DISCUSSION AND FUTURE WORK 

We presented two algorithms for reconstructing the intrin- 
sic correlations between distance-dependent quantities in 
apparent magnitude limited photometric redshift datasets. 
One was a generalization of the non-parametric Vmax 
method (Section [Sj, and the other used a maximum- 
likelihood approach (Section 

Both our reconstruction methods assume that the dis- 
tribution of photo-2 errors is known accurately. In practice, 
this means that spectroscopic redshifts are available for a 
subset of the data. The question then arises as to whether 
or not the number of spectra which must be taken to specify 
the error distribution reliably is sufficient to also provide a 
reliable (spectroscopic!) estimate of these scaling relations. 
If so, what is the basis for deciding that it is worth recon- 
structing these relations from the photo-z data? This is the 
subject of work in progress, although the methods presented 
in this paper assume that such reconstructions will indeed 
be necessary. E.g., if the spectra are not simply a random 
subset of the magnitude limited photometric sample, then it 
may be difficult to quantify and so correct for the selection 
effects associated with the spectroscopic subset. 

We used the size-luminosity relation in a mock cata- 
log which had realistic choices for the correlation to illus- 
trate the biases which are present and must be corrected 
if photometric redshift datasets are to provide reliable es- 
timates of galaxy scaling relations (Figures [2] and |3]). We 
showed that our iterative deconvolution scheme provides a 
simple and reliable correction of this bias (Figures [3H6]). 
Note that although we have illustrated our methods using 
a 2-dimensional distribution, the extension to n-correlated 
variables is trivial. 

Because our algorithm permits the accurate measure- 
ment of many scaling relations for which spectra were pre- 
viously thought to be necessary (e.g. the color-magnitude 
relation, the size-surface brightness relation, the Photomet- 
ric Fundamental Plane), we hope that our work will permit 
photometric redshift surveys to provide more stringent con- 
straints on galaxy formation models at a fraction of the cost 
of spectroscopic surveys. 

Our results may have other applications. For example, 
Bernard! (2007) has highlighted a bias associated with the 
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